R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

#Install packages

install.packages("BiocManager")
BiocManager::install("dada2", version = "3.11")

install.packages("qiime2R")
install.packages("tidyverse")
install.packages("biomformat")
install.packages('vegan')
install.packages('readr')
 

#Load packages

library("tidyverse")

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------------------------------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.2     v purrr   0.3.4
v tibble  3.0.3     v dplyr   1.0.2
v tidyr   1.1.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------------------------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library("dada2")

Loading required package: Rcpp

library("biomformat")
library('vegan')

Loading required package: permute
Loading required package: lattice
This is vegan 2.5-6

library('readr')
library("dplyr")     # To manipulate dataframes
library("readxl")    # To read Excel files into R
library("ggplot2")   # for high quality graphics
library("BiocManager")

Bioconductor version 3.11 (BiocManager 1.30.10), ?BiocManager::install for help
Bioconductor version '3.11' is out-of-date; the current release version '3.12' is available with R version '4.0'; see
  https://bioconductor.org/install

library("phyloseq")   
library("tidyr")
library("data.table")

data.table 1.13.0 using 2 threads (see ?getDTthreads).  Latest news: r-datatable.com

Attaching package: 㤼㸱data.table㤼㸲

The following objects are masked from 㤼㸱package:dplyr㤼㸲:

    between, first, last

The following object is masked from 㤼㸱package:purrr㤼㸲:

    transpose

#Merging tables for species abundance 1
#http://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/



###read in SILVA biom-taxonomy for SILVA ####################################################

taxonomy45 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/V4V5_SILVA_taxonomy.csv")
taxonomy68 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/taxonomy_SILVA_V6V8.csv")

head(taxonomy45)

#make a data.table
taxonomy45<-as.data.table(taxonomy45)
taxonomy68<-as.data.table(taxonomy68)
taxonomy45

taxonomy68

###read in otu table
otu45<- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/V4V5__ASV_table.csv", header = T) #,row.names=1
otu68<- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/V6V8__ASV_table.csv", header = T)
head(otu45) 

head(otu68)

#Convert to long format
otu_long45<-gather(otu45, Sample, Counts, "X122":"X83") #dashes got changed to dots for some reason
otu_long68<-gather(otu68, Sample, Counts, "X122.o":"X83.o") #dashes got changed to dots for some reason

#make a data. Table
otu_long45<-as.data.table(otu_long45)
otu_long68<-as.data.table(otu_long68)
otu_long45

otu_long68

###read in metadata file
metaD45<-read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/metadata_16S_V4V5_jonesSound2019.txt", sep = "")
metaD45<-as.data.table(metaD45)
metaD45$Sample<-as.character(metaD45$Sample)
metaD45

metaD68<-read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/metadata_16S_V6V8_jonesSound2019.txt", sep = "")
metaD68<-as.data.table(metaD68)
metaD68$Sample<-as.character(metaD68$Sample)
metaD68

# first set the keys to the common column:
setkey(otu_long68,Sample)
setkey(metaD68,Sample) 

setkey(otu_long45,Sample)
setkey(metaD45,Sample) 

# join the tables
total1_45 <- merge(otu_long45,metaD45,by="Sample")
total1_45

total1_68 <- merge(otu_long68,metaD68,by="Sample")
total1_68

#set key of total1_45 as the OTUID to join the taxonomy table
setkey(total1_45,OTUID)
setkey(taxonomy45,OTUID)

setkey(total1_68,OTUID)
setkey(taxonomy68,OTUID)

#merge to big table
total2_45 <- merge(total1_45,taxonomy45)
setkey(total2_45, Sample)
total2_45

total2_68 <- merge(total1_68,taxonomy68)
setkey(total2_68, Sample)
total2_68

#Final Joined table of SILVA
all_joined_SILVA45 <- metaD45[total2_45]
(all_joined_SILVA45)

all_joined_SILVA68 <- metaD68[total2_68]
(all_joined_SILVA68)

# Join table for Greengenes ####################################

#made above
otu_long45

otu_long68

metaD45

metaD68

total1_45

total1_68

taxonomyGG45 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/taxonomy_greengenes_V4V5.csv")
taxonomyGG68 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/taxonomy_greengenes_V6V8.csv")
taxonomyGG45 <- as.data.table(taxonomyGG45)
taxonomyGG68 <- as.data.table(taxonomyGG68)

setkey(taxonomyGG68,OTUID)
setkey(taxonomyGG45,OTUID)

total2_GG45 <- merge(total1_45,taxonomyGG45)
setkey(total2_GG45, Sample)
total2_GG45

total2_GG68 <- merge(total1_68,taxonomyGG68)
setkey(total2_GG68, Sample)
total2_GG68

#Final Joined table of GreenGenes
all_joined_GG45 <- metaD45[total2_GG45]
(all_joined_GG45)

all_joined_GG68 <- metaD68[total2_GG68]
(all_joined_GG68)

# Join table for RDP ####################################

#made above
otu_long45

otu_long68

metaD45

metaD68

total1_45

total1_68

taxonomyRDP45 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/taxonomy_RDP_V4V5.csv")
taxonomyRDP68 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/taxonomy_RDP_V6V8.csv")
taxonomyRDP45 <- as.data.table(taxonomyRDP45)
taxonomyRDP68 <- as.data.table(taxonomyRDP68)

setkey(taxonomyRDP68,OTUID)
setkey(taxonomyRDP45,OTUID)

total2_RDP45 <- merge(total1_45,taxonomyRDP45)
setkey(total2_RDP45, Sample)
total2_RDP45

total2_RDP68 <- merge(total1_68,taxonomyRDP68)
setkey(total2_RDP68, Sample)
total2_RDP68

#Final Joined table of GreenGenes
all_joined_RDP45 <- metaD45[total2_RDP45]
(all_joined_RDP45)

all_joined_RDP68 <- metaD68[total2_RDP68]
(all_joined_RDP68)

NA
NA
NA
NA

#####Proportion plot for each sample with both primer, based on phylum (?) ######
p45 <- ggplot(all_joined_SILVA45, aes(x = Counts, y = Sample, fill = Phylum)) + geom_bar(position = "fill", stat = "identity")

Error in ggplot(all_joined_SILVA45, aes(x = Counts, y = Sample, fill = Phylum)) : 
  could not find function "ggplot"

  
######TABLE 1: table of prevalence of bacterial counts according to each primer#####
#bacteria on x, count on y, two different graphs or two coloured bars on each graph

otu_tax_45 <- merge(otu_long45,taxonomy45,by="OTUID")
otu_tax_45

otu_tax_68 <- merge(otu_long68,taxonomy68,by="OTUID")
otu_tax_68

# Basic barplot
p<-ggplot(data = otu_tax_45, aes(y=Phylum, x=Counts/1000)) +
  geom_bar(stat="identity",position = position_dodge(width=1.5))
p

plot2 <- ggplot(NULL, aes(y=Phylum, x=Counts/1000)) + 
      geom_bar(data = otu_tax_45, stat = "identity", colour="red",width=0.4, position = position_dodge(width=0.5)) +
      geom_bar(data = otu_tax_68, stat = "identity", colour="blue",width=0.4, position = position_dodge(width=0.5))
plot2

all_joined_SILVA45

all_joined_SILVA68

### V4V5 ###
orderCount45<- select(all_joined_SILVA45, Order, Counts) #count by Order.
orderCount45<- orderCount45[order(orderCount45$Counts)]

chloroTable45 <- filter(orderCount45, Order == "Chloroplast") #Select only rows that are chloroplasts
chloroSum45 <- sum(chloroTable45$Counts)
chloroSum45 #this is the number of sequences associated with chloroplasts by V4V5

[1] 78158

### V6V8 ###
orderCount68<- select(all_joined_SILVA68, Order, Counts) #count by Order.
orderCount68<- orderCount68[order(orderCount68$Counts)]

chloroTable68 <- filter(orderCount68, Order == "Chloroplast") #Select only rows that are chloroplasts
chloroSum68 <- sum(chloroTable68$Counts)
chloroSum68 #this is the number of sequences associated with chloroplasts by V4V5

[1] 57941

#TODO
#Counts of archaea members SILVA

#Database like commands for tables!!!
#https://stackoverflow.com/questions/12353820/sort-rows-in-data-table-in-decreasing-order-on-string-key-order-x-v-gives-er

#For V4V5
all_joined_SILVA45
sortTable45 <- all_joined_SILVA45[order(rank(Kingdom),Phylum)] #sort by kingdom name
Archaeas45 <- filter(sortTable45, Kingdom == "Archaea") #Select only rows where Kingdom is Archaea.
Archaeas45

#Simple bar graph
archPerSample45 <-ggplot(Archaeas, aes(x=Counts, y=Sample)) +
  geom_bar(stat="identity") +
  theme_light()
archPerSample45


#For V6V8
all_joined_SILVA68
sortTable68 <- all_joined_SILVA68[order(rank(Kingdom),Phylum)] #sort by kingdom name
Archaeas68 <- filter(sortTable68, Kingdom == "Archaea") #Select only rows where Kingdom is Archaea.
Archaeas68

#Simple bar graph
archPerSample68 <-ggplot(Archaeas68, aes(x=Counts, y=Sample)) +
  geom_bar(stat="identity") +
  theme_light()
archPerSample68

orderNames_S68 <- c("Nitrosopumilales")
orderSums_S68 <- c() #sums of each of the orders

#loop to find the sums of each order and add to the vector
for (oName in orderNames_S68) {
  curOrderTable <- filter(SILVA_table68, Order == oName) #Select only rows with the current order name from the full table
  curOrderTable <- select(curOrderTable, Order, Counts) #makes a sub table of only the current table's order name and counts
  
  curOrderSum <- sum(curOrderTable$Counts) #sums all the counts of that order
  orderSums_S68 <- c(orderSums_S68, curOrderSum) 
}
(orderSums_S68)

[1] 8226

orderDataFrame_S68 <- data.frame(orderNames_S68, orderSums_S68)
orderDataFrame_S68

topArchGraph_S68 <-ggplot(orderDataFrame_S68, aes(x=orderSums_S68, y=orderNames_S68)) +
  geom_bar(stat="identity") +
  theme_light()
topArchGraph_S68

#TODO
#Phyloseq to get alpha diversity
library(phyloseq)

otu_mat<- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/V4V5__ASV_table.csv", header = T, row.names = 1) #,row.names=1
tax_mat <-  read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/V4V5_SILVA_taxonomy.csv", header = T, row.names = 1)
samples_df<-read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/metadata_16S_V4V5_jonesSound2019.txt", sep = "")

head(tax_mat)
head(otu_mat)

#assign row names
row.names(otu_mat) <- otu_mat$OTUID
#row.names(tax_mat) <- tax_mat$OTUID
#row.names(samples_df) <- samples_df$Site

#Transform into matrixes otu and tax tables (sample table can be left as data frame)
otu_mat <- as.matrix(otu_mat)
#tax_mat <- as.matrix(tax_mat)

class(tax_mat)

#Transform to phyloseq objects
class(otu_mat) <- "numeric" #NEED to convert to numeric in this way. Have to run whole chunk.

OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)

head(OTU)
head(TAX)

#taxa_names(OTU)
sample_names(OTU)
sample_names(TAX)

taxa_names(OTU)
taxa_names(TAX)

#Convert phyloseq
Perma <- phyloseq(OTU, TAX, samples)
Perma

#Making Bubble plots (TODO)

#V4V5 -----------------------------------------------------------------------------------------------------------------

###read in otu table of V4V5
otu45<- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/V4V5__ASV_table.csv") #,row.names=1
head(otu45) 
#Convert to long format
otu_long45<-gather(otu45, Sample, Counts, "X122":"X83") #dashes got changed to dots for some reason
#make a data. Table
otu_long45<-as.data.table(otu_long45)
otu_long45

###read in V4V5_SILVA_taxonomy.csv
taxonomy45 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/V4V5_SILVA_taxonomy.csv")
head(taxonomy45)
#make a data.table
taxonomy45<-as.data.table(taxonomy45)
taxonomy45

###read in metadata file
metaD45<-read.csv2("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/metadata_16S_V4V5_jonesSound2019.txt")
metaD45<-as.data.table(metaD45)
head(metaD45)


write.csv(otu_long45, "G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V4V5/for_plot45.csv")
#remove 0s
for_plot45 <- filter(otu_long45, Counts!= 0)
#for relative abundance:
#change "Name" for the name of your column with group description
for_plot45$type<-as.factor(for_plot45$type)
#order
for_plot45$type<-factor(for_plot45$type, levels=c( "river", "natural biofilm", "artificial biofilm"))
#Now order your x axis (site names) <-- edit accordingly.
for_plot45$Site <- factor(for_plot45$Site, levels = c( "Sif_1","SI_2","SI_3", "S4_1", "Big_1", "BH_2", "BH_3", "Ram_1","RAM_3", "Bapt_1", "BP_2", "BP_3", "Nord_1", "N_2", "N_3", "Rose_1", "RMS_1", "RMS_2", "RM_3", "MM_1", "Mod_1", "MM_2", "MM_3", "TM_1", "T_2", "SM_2", "White_1", "WM_2", "WM_3", "GM_1", "Sturg_1", "ST_2", "Verm_1", "VM_1", "VA_1", "VM_2", "VM_3", "V3_1"))

#sweet- now plot (note that this is for raw counts; change for relative abundance)
theme_set(theme_bw())
col=c( "green","orange", "purple")## <- change accordingly, if desired-- see http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
g <- ggplot(for_plot45, aes(x=Sample, y=Sample)) 
g + geom_point(aes(size=Counts)) +
  scale_colour_manual(values=col)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=9))+
  theme(axis.text.y = element_text( size=9))+
  theme(axis.title.x=element_text(size=12))+
  theme(axis.title.y=element_text(size=12))
#export both bubble plots as PDF.


#V6V8 -----------------------------------------------------------------------------------------------------------------

###read in otu table of V6V8 
otu68<- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/V6V8__ASV_table.csv") #,row.names=1
head(otu68) 
#Convert to long format
otu_long68<-gather(otu68, Sample, Counts, "X122.o":"X83.o") #dashes got changed to dots for some reason
#make a data. Table
otu_long68<-as.data.table(otu_long68)
otu_long68

###read in biom-taxonomy
taxonomy68 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/taxonomy_SILVA_V6V8.csv")
head(taxonomy68)
#make a data.table
taxonomy<-as.data.table(taxonomy68)
taxonomy68

###read in metadata file
metaD68<-read.csv2("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/metadata_16S_V6V8_jonesSound2019C.csv") 
metaD68<-as.data.table(metaD68)
head(metaD68)



write.csv(otu_long68,"G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/for_plot68.csv")
v68 <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Real_Content/Data_Files/V6V8/for_plot68.csv")
v68

#remove 0s
for_plot68 <- filter(otu_long68, Counts!= 0)
#for relative abundance:
#change "Name" for the name of your column with group description
for_plot68$type<-as.factor(for_plot68$type)
#order
for_plot68$type<-factor(for_plot68$type, levels=c( "river", "natural biofilm", "artificial biofilm"))
#Now order your x axis (site names) <-- edit accordingly.
for_plot68$Site <- factor(for_plot68$Site, levels = c( "Sif_1","SI_2","SI_3", "S4_1", "Big_1", "BH_2", "BH_3", "Ram_1","RAM_3", "Bapt_1", "BP_2", "BP_3", "Nord_1", "N_2", "N_3", "Rose_1", "RMS_1", "RMS_2", "RM_3", "MM_1", "Mod_1", "MM_2", "MM_3", "TM_1", "T_2", "SM_2", "White_1", "WM_2", "WM_3", "GM_1", "Sturg_1", "ST_2", "Verm_1", "VM_1", "VA_1", "VM_2", "VM_3", "V3_1"))

#sweet- now plot (note that this is for raw counts; change for relative abundance)
theme_set(theme_bw())
col=c( "green","orange", "purple")## <- change accordingly, if desired-- see http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
g <- ggplot(for_plot68, aes(x=Sample, y=Sample))#colour=type
g + geom_point(aes(size=Counts)) +
  scale_colour_manual(values=col)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size=9))+
  theme(axis.text.y = element_text( size=9))+
  theme(axis.title.x=element_text(size=12))+
  theme(axis.title.y=element_text(size=12))
#export both bubble plots as PDF.

#Data processing with Pyloseq:https://vaulot.github.io/tutorials/Phyloseq_tutorial.html


otu_mat <-read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Example_Content/all_OTU_NMDSnu.csv", row.names=1) 
tax_mat <- read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Example_Content/biom-taxonomy.csv", row.names=1)
samples_df  <- as.data.frame(read.csv("G:/Desktop/Folders/Edu/2] EAS/Thesis/Example_Content/all_OTU_NMDS.csv"))

head(otu_mat)
head(samples_df)

#assign row names
row.names(otu_mat) <- otu_mat$OTUID
row.names(tax_mat) <- tax_mat$OTUID
row.names(samples_df) <- samples_df$Site

#Transform into matrixes otu and tax tables (sample table can be left as data frame)
otu_mat <- as.matrix(otu_mat)
tax_mat <- as.matrix(tax_mat)

class(otu_mat)

#Transform to phyloseq objects
class(otu_mat) <- "numeric" #NEED to convert to numeric in this way. Have to run whole chunk.

OTU = otu_table(otu_mat, taxa_are_rows = TRUE)
TAX = tax_table(a)
samples = sample_data(samples_df)

nrow(OTU)
nrow(TAX)
nrow(samples)

#Convert phyloseq
Perma <- phyloseq(OTU, TAX, samples)
Perma